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The recent solution to the fermion sign problem allows, for the first time, the 
use of cluster algorithm techniques to compute certain fermionic path integrals. 
To illustrate the underlying ideas behind the progress, a cluster algorithm is 
. constructed to study the chiral phase transition in a strongly interacting stag- 

gered fermion model with an arbitrary mass term in 3 + 1 dimensions. Unlike 
^ ■ conventional methods there is no difficulty in the cluster method to approach 

. the chiral (massless) limit. Results using the new algorithm confirm that the 

I chiral transition falls under the expected universality class. 
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I. Introduction 

The past few years has been a time of remarkable progress in our understanding 
O"! of lattice fermions and their abihty to reproduce the chiral properties of the continuum 



^ ' theory It is now possible to formulate vector-like gauge theories on the lattice with 

all its chiral symmetry intact Our understanding of non-perturbative formulations 
^ ■ of chiral gauge theories is also advancing at a remarkable rate as is reflected by many 

. recent contributions 0. On the other hand, our ability to perform reliable (numerical) 

calculations in many theories involving strongly interacting fermions and especially in the 
chiral limit has remained primitive. In particular algorithms based on hybrid molecular 
dynamics, that are popular for gauge theories like QCD, suffer either from sign problems 
or critical slowing down in the chiral limit. Unless better numerical methods are found, 
it is likely that we will be forced to perform most of our calculations far from the chiral 
limit and hence not gain from the progress in our understanding of chiral symmetry. 

Numerical methods for fermionic systems are especially difficult due to the Pauli- 
exclusion principle. The Boltzmann weight of the statistical mechanics problem describing 
such systems can often become negative thereby invalidating most Monte-Carlo methods. 
This is usually referred to as the sign problem. In other cases one is often forced to 
work with bosonic variables with non-local Boltzmann weights which arise through the 
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fermion determinant. Such expressions are not easy to handle and there is no known way 
to overcome problems like critical slowing. 

Recently a new numerical method to perform fermionic path integrals has emerged 
and is based on a rather novel form of bosonization. Instead of the Grassmann variables, 
the method starts by representing fermions through their occupation numbers. Such 
representations, which have been studied in the past were not considered useful 

due to sign problems. Over the past year it has become clear that at least in a class of 
models it is possible to design cluster algorithms which automatically solve the fermion 
sign problem and provide an alternative way to compute fermionic path integrals. In a 
sense the algorithm suggests a new way to bosonize the model in terms of dynamics of 
clusters. The new method works in both relativistic ^ and non-relativisticP] models. 
Since cluster algorithms are known to beat critical slowing down in a variety of bosonic 
models]^, |lOl, the new algorithms have the ability to work efficiently even in the presence 
of long correlation lengths. In particular, they encounter no difficulties in the limit where 
the fermions become massless. 

The present article illustrates the main features of the new approach through an 
example of a four Fermi model involving staggered fermions. This model undergoes a 
Zi chiral transition from a high temperature chirally symmetric phase to the broken 
phase at low temperatures and was studied originally in the chiral limit in |]^. Here 
the earlier work is extended to include a fermion mass. This extension turns out to be 
quite useful because on a large but finite lattice all the interesting properties of the chiral 
condensate, as an order parameter of the transition, emerge only when a tiny mass can 
be added to the system. In the next section the partition function of the model is written 
in the occupation number basis and the origin of the fermion sign problem is reviewed. 
In section 3 the partition function is rewritten in terms of connected fermion world-line 
configurations that arise due to the introduction of local "bond" variables. In the new 
variables the solution to the sign problem emerges naturally and leads to a very elegant 
cluster algorithm. In section 4, the algorithm and some numerical results are discussed. 
Section 5 contains conclusions along with a description of on-going work and points to 
new directions for the future. 



II. The Model 



The model considered in this article involves interacting staggered fermions hopping 
on a 3-d cubic spatial lattice with V = sites x {L even) and with anti-periodic spatial 
boundary conditions. The dynamics of the fermions is described through the Hamilton 
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operator 

H = J2hx,i + mJ2Px (1) 

x,i X 

where the term 

hx, = IvxA'^t^x+l + + (^^^. - \)iK+i^x+l - I), (2) 

couples the fermion operators at the lattice sites x and x + i, where i is a unit- vector in 
the z-direction and the mass term 

p, = (-ir+^^+^«(vl/+v|/,-l/2), (3) 

is a single site operator. The fermion creation and annihilation operators \E'^ and used 
in the above equations satisfy the standard anti-commutation relations 

{VI/+ VI/+} = {M/,, VI/J = 0, {VI/+, ^y} = 6xy. (4) 

Further rj^^i = 1, ?7x,2 = and rj^^s = (— g^j^g ^j^g standard staggered fermion 

sign factors. In the chiral limit (m = 0), the above Hamiltonian was discussed in 0. 

In addition to the usual U{1) fermion number symmetry, there are a number of dis- 
crete symmetries of the staggered fermion Hamiltonian as discussed in |^ . For example, 



at m = the Hamiltonian is invariant under shifts of \E'|jr and by one lattice spacing 
in the x^ direction. The mass term breaks this Z2 symmetry which is known to be a 
subgroup of the well known chiral symmetry of relativistic massless fermions. Here this 
symmetry is broken spontaneously at low temperatures while thermal fluctuations restore 
it at high temperatures. The associated second order phase transition belongs to the 
universality class of the 3d ising model as discussed in The extension discussed here, 
allows one to add a tiny mass and hence study the transition using the chiral condensate. 

In order to construct the path integral representation of the partition function the 
Hamilton operator is first decomposed into seven terms H = Hi + H2 + ■■■ + Hq + with 

Hi= ^»+3 = H for i = 1,2,3; = mJ^Px- (5) 

x = (xi,X2,x^) x = (xi,X2,X3) ^ 

XiSven Xiodd 

This decomposition is exactly the same as the one used in for the m = case, where 
Hj was absent. The various steps in the construction of the partition function from here 
on is exactly the same as in 0. Leaving the details to that paper only the essential points 
are sketched here. The Suzuki- Trotter formula leads to the expression 

Zf = Tr[exp(-/3if)] = hm Tr[e-^^^e-^^^/^e-^^^e-^^^/^..e-^^''e-^^^/^]^^ (6) 
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for the fermionic partition function at inverse temperature (3, where e = /3/M is the 
lattice spacing in the Euchdean time direction. Spreading Hj symmetrically with each 
i^i, i = 1, ..6 is not necessary except that it adds some symmetry to the algorithm. Further, 
in the arguments given below, M is assumed to be finite and that the true partition 
function is obtained from an extrapolation of the data obtained from a series of simulations 
at larger and larger M. Fortunately it is also possible to formulate the cluster algorithm 
in the time continuum limit as demonstrated in 



t=2£ 



t=£ 



t=0 



► 

X 

FIG. 1. A fermion world- line configuration with Signf[n] = — 1 on a one dimensional periodic 
spatial lattice of size L = 8 where the time slice t = 3e is the same as t = 0. 

It is well known that in the occupation number basis a fermionic partition function 
can be written as a sum over fermion world-line configurations[^]. Since the fermionic 
creation operators anti-commute, the information that a subset of the lattice sites are 
occupied defines the state of the system only up to a sign. A complete definition of 
the occupation number basis states requires an ordering of the entire lattice which then 
specifies the order in which the fermions on the lattice are created. Taking this ordering 
into account it is possible to write the partition function as a sum over configurations of 
fermion occupation numbers n{x, t) = 0, 1 on a (3 -|- l)-d space-time lattice of points {x, t) 
where t = en/12, n = 0, 1,2, ...(12M — 1) labels the time slices. The trace imposes the 
periodicity constraint n{x,0) = n{x,P). A typical configuration in one spatial dimension 
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is shown in figure 1, where the shaded regions depict either a two site interaction due to 
exp(— /la; j) or a single site interaction due to exp(— mp^;). 

The Bohzmann weight of each configuration n is a product of transfer matrix ele- 
ments associated with the configuration of fermions on each shaded region. These elements 
are shown in figure 2. Representing the magnitude of the Boltzmann weight by Win] and 

Non-zero elements due to exp(— /la- j) Non-zero elements due to exp^—p^) 




FIG. 2. The figures (a) through (j) illustrate the non-zero transfer matrix elements. Weight 
of (a),(b) is exp(— e/4); weight of (c),(d) is exp(e/4) cosh(e/2); weight of (e),(f) 
is S exp(e/4) sinh(e/2); weight of (g),(j) is exp(— me/12) and weight of (h),(j) is 
exp(me/12). The factor S is a product of local sign factors 7]x^i and non-local sign 
factors that arise due to anti-commutation relations. Further, xi + X2 + X3 = even(odd) 
defines an even(odd) site. 

the sign, which takes into account the various anti-commutation relations, by Sign[n] the 
partition function takes the form 

Zf = Y.Sign[n]W[n]. (7) 

n 

Clearly Sign[n] is the product of the sign factors E that appear in the transfer matrix 
elements associated to fermion hops represented through figures 2(e) and 2(f). When 
the fermion hops from x to x + i ot vice versa due to the action of exp(— /i^. j) it picks 
up a product of local sign factors due to terms like rij;^i as well as a string of non-local 
signs that arise due to anti-commutation relations involved when the fermion has to cross 
other fermions on the ordered lattice while reaching its destination. The evaluation of 
this non-local part of S is rather tedious. Separating the two parts, one can define 
Sign[n] = Sign£[r;,]Signjj[r;,], where Signj[r;,] comes from the product of non-local parts of S 
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and is purely fermionic, where as Sign{,[n] is the one that comes from the local parts and 
may arise even in bosonic models. Fortunately Sign£[n] has a topological meaning. The 
occupied lattice sites define fermion world-lines which are closed around the Euclidean 
time direction. However, during their Euclidean time evolution fermions can interchange 
their positions, and the fermion world-lines thus define a permutation of particles. The 
Pauli exclusion principle dictates that the Signj[n] is just the sign of that permutation. 
Further, Sign(5[n] receives an extra minus-sign for every fermion that hops across a spatial 
boundary due to anti-periodic boundary conditions. 




FIG. 3. The figure on the left shows a connected fermion world-line configuration in one space 
and one time dimension. The figure on the right is obtained after flipping a cluster 
(which is shown) from the configuration on the left. The resulting new connected 
fermion world-line configuration is the "reference" configuration of the model. 

III. Bond Variables and Meron Clusters 

When the partition function is written as a sum over fermion world-line configura- 
tions, the Boltzmann weight of each configuration is not guaranteed to be positive definite. 
This leads to the well known fermion sign problem. It was discovered recently that some- 
times it is possible to solve this problem when the Boltzmann weight of each fermion 
world line configuration is written as a sum of Boltzmann weights describing new type 
of configurations, referred to as "connected" fermion world line configurations, obtained 
by the introduction of bond variables that establish connections with in the lattice. An 
example of this new type of configuration is shown in figure 3. In a class of models. 
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FIG. 4. The figure shows how to rewrite the weight of each configuration of figure |2|, in terms 
of new configurations that have bonds in addition to fermions. 



7 



including the one considered in this article, it is possible to regroup these new type of 
configurations such that the sum of the weights of each group is positive or zero, thus 
canceling all negative weights. The regrouping is accomplished by identifying all possible 
connected fermion world-line configurations obtained by fiipping clusters of connected 
lattice sites. By fiipping a lattice site one means that if the site is filled (with a fermion) 
it is emptied and vice versa. The result of a cluster fiip is also shown in Figure 3. 

In order to calculate the Boltzmann weight of every new type of configurations each 
interaction plaquette of in the old configuration given in figure 2 is rewritten as a sum 
of new configurations that involve bonds in addition to the fermions. The sign factors S 
are ignored at this step and will be considered during the regrouping step. Further, the 
bonds are introduced so that they obey the property of a flip symmetry. This means that 
the magnitude of the Boltzmann weight of connected fermions world line configurations 
do not change when connected sites are fiipped. This property will be useful to identify 
configurations with equal weight but opposite signs during the regrouping step. For the 
two site interactions the weights of the new bond variables is shown at the top of figure 4. 
For the single site interaction due to the mass terms one has to distinguish between even 
and the odd sites since the weights differ on these two sites. However, now by introducing 
two types of bonds the property of flip symmetry can be maintained as shown in the 
bottom part of flgure 4. The difference between the two bonds is that when the sites 
connected by the "dashed" bond are flipped, an extra sign change has to be taken into 
account. These signs must be multiplied to the Sign[n] and the overall sign should be 
taken into account at the regrouping stage. 

The above procedure gives a new expression for the path integral in terms of con- 
nected fermion world-line conflgurations. 

Zf = Y.Sign[C]W[C] (8) 
c 

Each conflguration C now deflnes a set of clusters in addition to the fermion occupation 
numbers as shown in flgure 0. For the present model the clusters form closed loops. The 
weight W[C] of each conflguration can be read off using the rules of flgure ^. The Sign[C] 
is the product of Sign[n] for the old conflguration obtained by ignoring the bond and 
the extra negative sign factors that may arise from fllled dashed-bonds on even sites and 
empty dashed-bonds on odd sites. A great simpliflcation occurs when the flip symmetry 
imposed during the bond introduction is used. If a conflguration C has Nc clusters, then 
all the 2^'^ conflgurations, obtained by flipping the clusters independently, have the same 
weight iy[C]. This degeneracy can be used to regroup the conflgurations and write 

Zf = Y.Si^]W[C] (9) 

c 
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where Sign[C] is the average of Sign[C] over these 2^^ configurations. The usefulness of 
eq. can be appreciated when Sign[C] obeys the following two properties: 

(1) Sign[C] = Hi^i Sign[cj], where Sign[cj] is the sign of the ith cluster configuration 
labeled q. 

(2) For every configuration the individual clusters can be fiipped to a reference config- 
uration cf^ such that Sign[c^'^^] = 1. 



It is then easy to show that Sign[C] = 0, if there is at least a single cluster whose fiip 
changes the sign of the configuration, i.e., there is a cluster q such that Sign[c"°'^~''^^] = —1. 
Such clusters are referred to as merons. A configuration with no merons will then have 
Sign[C] = 1. There are models in which Sign[C] has the above two properties and the 
present model falls in this class. Although the first property of Sign[C] is difficult to 
establish, it is easy to show that every cluster in the present model can always be fiipped 
such that all even sites are empty and odd sites are filled in the present model. This leads 
to the staggered configuration shown on the right in figure |^. 

The new expression for the partition function derived in eq. shows that fermion 
dynamics in the present model is equivalent to a statistical mechanics of clusters. Each 
cluster has two degrees of freedom related to the fermions. The Pauli-exclusion principle 
is encoded in the fact that meron clusters contribute zero weight to the partition function. 
In order to determine if a cluster is a meron or not, one has to understand its topology. 
If Tih is the number of hops of the cluster to the neighboring lattice site on the same time 
slice, is the number of local signs that the cluster encounters during the hops, ra^ is the 
temporal winding and is the number of dashed bonds in the cluster, then the cluster 
is a meron if (n^ + nh/2 + + n^) is an even integer. Cluster algorithms based on the 
concept of merons was originally invented in and are referred to as the meron cluster 
algorithms. 



IV. Algorithm and Results 



It is easy to construct cluster algorithms that generate clusters with weight W[C] 
obtained using the rules given in figure ^ This makes it possible to simulate the modified 
model obtained by ignoring the sign factors. The partition function of the resulting model, 
denoted by Zb, can be expressed as a sum over contributions from various meron sectors, 
i.e., Zb = YjN^n where N denotes the number of merons in a cluster configuration. On 
the other hand the actual partition function Zj = Zo. Since a typical configuration of 
Zb is filled with merons, it is exponentially difficult to find zero meron sectors and the 
cluster algorithm is still inefficient for simulating the dynamics of Zf. A solution to the 
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FIG. 5. The figure shows the chiral condensate {ipip) as a function of the inverse temperature 
j3 at m = 0.001 for various spatial volumes. 



problem of avoiding configurations with meron clusters is necessary. Fortunately, using a 
local Metropolis accept-reject step, it is possible to introduce a re-weighting factor that 
suppresses higher meron sectors. For example the Metropolis step can modify the weight 
of cluster configurations to VF[C]/p^, where p > 1 and represents the number of merons 
in the configuration. This modifies Zb but leaves Zf unchanged! This Metropolis step 
then increases the efficiency of the algorithm by an exponential factor. 

A variety of quantities can be calculated with the new algorithm. The first step, 
however, is to find expressions for these quantities in terms of cluster variables. Operators 
that are diagonal in the occupation number basis are easy to calculate, although more 
general operators can also be evaluated. For some quantities it is possible to find analytic 
expressions for the average over all the cluster fiips, which are referred to as improved 
estimators. The chiral condensate in the present problem, which is also the order param- 
eter for the chiral phase transition, is one such quantity. In the operator language the 
condensate is defined as 

where O — J2xPx- The improved estimator on the other hand is given by 

. - , 1 {p Sizc(Cmeron) Sn,i) , v 

2 {dN,o} V jj 
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where is the number of merons in the cluster configuration generated by the algorithm 
and Size(cmeron) is the size of the meron cluster. The extra re-weighting factor p takes 
into account the modifications due to the Metropolis step. It is easy to show that on a 
finite lattice in the m = limit there are no single meron sectors and hence the chiral 
condensate simply vanishes. On the other hand, as seen in figure an extremely tiny mass 
m = 0.001 shows convincingly the effects of spontaneous symmetry breaking. The critical 
value of = 0.948 obtained through finite size scaling analysis in appears consistent 
with the present results. In the chiral susceptibility was measured in the chiral limit 
its behavior near the transition was shown to be consistent with the universality class of 
the 3-d ising model. A more detailed analysis of the data from the present study will be 
presented elsewhere. 



V. Conclusions 



The fermion cluster algorithm constructed in this article is an example of a novel way 
to compute fermionic path integrals. The sign problems that arise in the new method are 
solved using the concept of a meron which makes it easy to match all the configurations 
with negative Boltzmann weight with configurations with an equal but positive Boltzmann 
weight. Some interesting features emerge with the new algorithm. For the first time it is 
possible to treat fermions and bosons on an equal footing in simulations fl^. Further it 
is possible to approach chiral limits more easily as discussed here. 

A number of applications of the new techniques are presently being studied. For 
example it appears possible to find models with continuous chiral symmetry in which 
the sign problem can be completely solved. Chiral symmetry breaking in such a system 
will lead to the existence of massless Goldstone particles in the chiral limit. Apart from 
clarifying issues of universality the presence of the low mass particles can help in studying 
resonance physics that become difficult with conventional algorithms. Applications to 
non-relativistic many fermion problems is another topic where the present techniques are 
directly applicable. It is also likely that solutions to Hubbard type models will emerge 
with these techniques. 

Adding gauge fields introduces new negative signs. It will be exciting to find solu- 
tions to this class of sign problems using cluster methods. There is evidence from strong 
coupling limits that a regrouping of configurations that differ in both gauge and fermion 
content can help solve such problems. This question is presently under investigation. 
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